
// register flip
#define rf(a, b)                              \
        do {                                            \
                float4 x = a;                           \
                int4 c = a> b;                          \
                a.s012 = select(b, a, c).s012;          \
                b.s012 = select(x, b, c).s012;          \
        } while (0)


__kernel void Median(__read_only image2d_t src, __write_only image2d_t dst){
    
    
        int gx = get_global_id(0), gy = get_global_id(1);
	
		__private sampler_t smp = CLK_NORMALIZED_COORDS_FALSE | CLK_ADDRESS_CLAMP_TO_EDGE | CLK_FILTER_NEAREST;
		
        float4 s0 = read_imagef(src, smp, (int2) { gx-1, gy-1 });
        float4 s1 = read_imagef(src, smp, (int2) { gx, gy-1 });
        float4 s2 = read_imagef(src, smp, (int2) { gx+1, gy-1 });
        float4 s3 = read_imagef(src, smp, (int2) { gx-1, gy });
        float4 s4 = read_imagef(src, smp, (int2) { gx, gy });
        float4 s5 = read_imagef(src, smp, (int2) { gx+1, gy });
        float4 s6 = read_imagef(src, smp, (int2) { gx-1, gy+1 });
        float4 s7 = read_imagef(src, smp, (int2) { gx, gy+1 });
        float4 s8 = read_imagef(src, smp, (int2) { gx+1, gy+1 });
		
		// credit for the branchless sort idea form google code project  "socles"
		// uses the median of medians algerithom 
		// on my machien it is 20% faster then with smart sorting statements
		
        //0
        rf(s0, s1);
        rf(s1, s2);
        rf(s0, s1);

        // 1
        rf(s3, s4);
        rf(s4, s5);
        rf(s3, s4);

        // 2
        rf(s6, s7);
        rf(s7, s8);
        rf(s6, s7);
		
        // final sort
        rf(s1, s4);
        rf(s4, s7);
        rf(s1, s4);
        

        write_imagef(dst, (int2) (gx, gy), s4);
}
